# for seekgene CRC human sample
library(Seurat)
library(tidyseurat)
source('00_util_scripts/mod_seurat.R')
library(tidyverse)

# seekgene starsolo velocyto ----------
skg_path <- '/home/supervisor/mist/gj_seekgene_crc/empty_intron/'

skg_raw <- list.files(skg_path, pattern = 'raw_full$', include.dirs = TRUE,
                      full.names = TRUE, recursive = TRUE) |>
  str_subset('Gene')

skg_final <- list.files(skg_path, pattern = 'filtered_full$', include.dirs = TRUE,
                        full.names = TRUE, recursive = TRUE) |>
  str_subset('Gene')

sobj_raw <- skg_raw |>
  set_names(c('crc221125','crc230304','crc230309')) |>
  Read10X() |>
  CreateSeuratObject(min.features = 1)

# do not filter any genes since soupx requires them
sobj_sg <- skg_final |>
  set_names(c('crc221125','crc230304','crc230309')) |>
  Read10X() |>
  CreateSeuratObject(min.features = 200)

# save original raw matrix for soupx correction
# ~2 min
system.time(raw_list <- sobj_raw |>
  SplitObject(split.by = 'orig.ident') |>
  map(LayerData, 'counts'))

# filter raw mtx with R function ---------
# ~7 min
celldrop <- raw_list |>
  map(DropletUtils::emptyDrops, .progress = TRUE) |>
  map(as_tibble, rownames = '.cell') |>
  map(filter, FDR < .01) |>
  map(pull, .cell) |>
  list_c()

# keep only non-empty drops
sobj_raw <- sobj_raw |>
  filter(nFeature_RNA >= 200 & .cell %in% celldrop)

sobj_raw$mitoRatio <- sobj_raw |> PercentageFeatureSet("^MT-")

sobj_raw <- sobj_raw |>
  filter(mitoRatio < 10)

sobj_raw <- sobj_raw |> quick_process_seurat()

cluster_list <- sobj_raw |>
  as_tibble() |>
  select(.cell, orig.ident, seurat_clusters) |>
  nest(data = -orig.ident) |>
  pull(data) |>
  map(\(x)set_names(x$seurat_clusters, x$.cell))

sobj_list <- sobj_raw |>
  SplitObject(split.by = 'orig.ident') |>
  map(LayerData, 'counts')

# or use starsolo filtered mtx ---------
sobj_sg$mito.ratio <- sobj_sg |> PercentageFeatureSet("^MT-")

sobj_sg |> ncol()
sobj_sg %<>%
  filter(mito.ratio < 10)
sobj_sg |> ncol()

sobj_sg %<>%
  quick_process_seurat()

cluster_list <- sobj_sg |>
  as_tibble() |>
  select(.cell, orig.ident, seurat_clusters) |>
  nest(data = -orig.ident) |>
  pull(data) |>
  map(\(x)set_names(x$seurat_clusters, x$.cell))

system.time(sobj_list <- sobj_sg |>
  SplitObject(split.by = 'orig.ident') |>
  map(LayerData, 'counts'))

## remove ambient RNA -----
des1 <- adjust_soup_solo(sobj_list[[1]], raw_list[[1]], cluster_list[[1]])

# cost ~2 min
system.time(desoup_list <- list(sobj_list, raw_list, cluster_list) |>
  pmap(adjust_soup_solo, .progress = TRUE))

desoup.mtx <- desoup_list[[1]] |>
  RowMergeSparseMatrices(desoup_list[2:3])

sobj <- desoup.mtx |>
  CreateSeuratObject(min.features = 200, min.cells = 3)

sobj$mito.ratio <- PercentageFeatureSet(sobj, "^MT-")

VlnPlot(sobj, 'mito.ratio')

sobj <- sobj %>%
  filter(mito.ratio < 10)

sobj <- quick_process_seurat(sobj)

sobj |> write_rds('CRC-I/data/seekgene.240923.rds')

## remove doublets ----------
# 0 task + progressbar MulticoreParam stuck > 10 min
# default SerialParam cost 3 min (result different from time to time)
# 3 worker Multicore cost 1.5 min 
system.time(
sobj.dbl <- sobj |>
  as.SingleCellExperiment() |>
  scDblFinder::scDblFinder(clusters = 'seurat_clusters',
                           samples = 'orig.ident',
                           BPPARAM = BiocParallel::MulticoreParam(3)) |>
  as.Seurat())
  
sobj.dbl$scDblFinder.class |> table()

sobj.dbl |>
  summarize(n = n(), .by = c(seurat_clusters, scDblFinder.class)) |>
  ggplot(aes(seurat_clusters, n, fill = scDblFinder.class)) +
  geom_col(position = "fill")

sobj.dbl |> DimPlot(cols = DiscretePalette(36))

sobj.dbl |> VlnPlot(c('mito.ratio','nCount_RNA'), group.by = 'seurat_clusters',
                    pt.size = 0)

sobj <- sobj.dbl |> filter(scDblFinder.class == 'singlet')

hpca <- celldex::HumanPrimaryCellAtlasData()

sobj <- mark_cell_type_singler(sobj, ref = hpca, new_label = 'hpca_main')

DimPlot(sobj, group.by = 'hpca_main',
        cols = DiscretePalette(36))

sobj <- sobj |>
  mutate(genotype = case_when(orig.ident == 'crc221125' ~ 'IT',
                              .default = 'TT')) |>
  filter(!str_detect(hpca_main, 'Endo|Epith|stem|muscle'))

write_rds(sobj, 'CRC-I/data/seekgene/crc2305.crEM.rds')

write_rds(sobj_raw, 'CRC-I/data/seekgene/crc2305.emtpydrop.rds')

# guo 2021 sided ---------
guo_sparse <- c('mtx','barc','feat') |>
  map(\(x)list.files('CRC-I/data/guo2021sided', pattern = x, full.names = TRUE)) |>
  pmap(ReadMtx, .progress = TRUE)

sobj_guo <- guo_sparse |>
  map2(c('left1','left2','left3','right1','right2','right3'), add_name_field) |>
  map(CreateSeuratObject, min.cells = 3, min.features = 200) |>
  purrr::reduce(merge)

sobj_guo %<>% JoinLayers()

sobj_guo$mito.ratio <- sobj_guo |>
  PercentageFeatureSet('^MT-')

sobj_guo |> VlnPlot('mito.ratio') 

sobj_guo <- sobj_guo |> filter(mito.ratio < 10)

sobj_guo <- sobj_guo |>
  quick_process_seurat()

## remove doublets ----------
# default SerialParam cost 2 min (result different from time to time)
# 6 worker Multicore cost 55s 
system.time(guo.scrbl <- sobj_guo |>
  mark_doublets())

system.time(res.dbl <- sobj_guo |>
              mark_doublets2())

sobj_guo %<>%
  left_join(res.dbl)

sobj_guo |>
  ggplot(aes(seurat_clusters, fill = scdbl.class)) +
  geom_bar(position = 'fill')

sobj_guo |> DimPlot(cols = DiscretePalette(36), label = T, label.box = T)

sobj_guo %<>% filter(scdbl.class == 'singlet')

hpca <- celldex::HumanPrimaryCellAtlasData()

sobj_guo <- mark_cell_type_singler(sobj_guo, hpca, new_label = 'hpca_main')

sobj_guo |> DimPlot(group.by = 'hpca_main',
                    cols = DiscretePalette(36))

immune_guo <- sobj_guo |>
  filter(str_detect(hpca_main, 'B|DC|Macro|Mono|Neutro|NK|T'))

immune_guo |> DimPlot(group.by = 'hpca_main',
                    cols = DiscretePalette(36))

immune_guo %<>%
  mutate(genotype = ifelse(str_detect(orig.ident, 'left'), 'IT', 'II'))

immune_guo |> dplyr::count(genotype)

immune_guo <- immune_guo |>
  quick_process_seurat(skip_norm = T)

immune_guo <- immune_guo |>
  mark_cell_type_singler(ref = hpca,
                         new_label = 'hpca_main')

immune_guo |> DimPlot(group.by = 'hpca_main') +
  ggtitle('Major tumor immune cell groups')

immune_guo |> write_rds('CRC-I/data/guo2021sided/guo2021_immune.rds')

# li2021liver ------------
mtx.liv <- Read10X('CRC-I/data/li_2021_liver')

sobj_liv <- mtx.liv |>
  CreateSeuratObject(names.field = 2, min.cells = 3, min.features = 200)

sobj_liv <- sobj_liv |> mutate(tissue = str_remove(.cell, '.+_'))

sobj_liv |> dplyr::count(tissue)

sobj_liv$mito.ratio <- sobj_liv |> PercentageFeatureSet("^MT-")

sobj_liv |> VlnPlot('mito.ratio', pt.size = 0)

sobj_liv <- filter(sobj_liv, mito.ratio < 10)

meta_liv <- "COL07\tIT\tGG
COL12\tII\tRR
COL15\tII\tGR
COL16\tIT\tGR
COL17\tIT\tGG
COL18\tII\tGG
" |> read_delim(col_names = c('orig.ident','FCGR2B.I232T','IGHG1.G396R'))

sobj_liv <- sobj_liv |> left_join(meta_liv)

# process 127k cells take ~11 min
sobj_liv <- sobj_liv |>
  quick_process_seurat(batch = c('orig.ident','tissue'))

sobj_liv |> DimPlot()

hpca <- celldex::HumanPrimaryCellAtlasData()

sobj_liv <- mark_cell_type_singler(sobj_liv,
                                   hpca,
                                   new_label = 'hpca_main')

sobj_liv |> DimPlot(group.by = 'hpca_main',
                    cols = DiscretePalette(10))

sobj_liv <- sobj_liv |> filter(!str_detect(hpca_main, 'Endo|Epi|stem'))

sobj_liv |> DotPlot('TNFRSF11B', group.by = 'hpca_main')

sobj_liv |>
  ggplot(aes(tissue, fill = hpca_main)) +
  geom_bar(position = 'fill')

sobj_liv |>
  filter(hpca_main == 'B_cell') |>
  FindMarkers(features = 'TNFRSF11B', group.by = 'tissue', ident.1 = 'PBMC',
              min.pct = 0)

sobj_liv |>
  filter(hpca_main == 'B_cell') |>
  VlnPlot(features = 'TNFRSF11B', group.by = 'IGHG1.G396R', pt.size = 0)

sobj_liv |>
  filter(str_detect(hpca_main, 'Epi')) |>
  FindMarkers(features = 'TNFRSF11B', group.by = 'seurat_clusters', ident.1 = 'GG',
              min.pct = 0, ident.2 = 'GR')

sobj_liv |>
  filter(str_detect(hpca_main, 'Epi')) |>
  VlnPlot(features = 'TNFRSF11B', group.by = 'IGHG1.G396R', pt.size = 0)

sobj_liv |>
  filter(str_detect(hpca_main, 'Epi')) |>
  bill.violin(features = 'TNFRSF11B', group.by = IGHG1.G396R)

sobj_liv |>
  as_tibble() |>
  select(.cell,'IGHG1.G396R')

sobj_liv |> write_rds('CRC-I/data/li_2021_liver/li2021tissue3.rds')

sobj_liv |>
  filter(tissue == 'CRC') |>
  write_rds('CRC-I/data/li_2021_liver/li2021crc.rds')

sobj_liv <- read_rds('CRC-I/data/li_2021_liver/li2021crc.rds')

sobj_liv |> DimPlot(cols = DiscretePalette(36))

sobj_liv %<>%
  quick_process_seurat(skip_norm = T)

# ~2 min: 1.3% doublet
# consider samples ~3.6 min: 0.7%
system.time(scrb.res <- sobj_liv |>
  mark_doublets())

scrb.res |>
  ggplot(aes(seurat_clusters, fill = scrublet_call)) +
  geom_bar(position = 'fill')

scrb.res |>
  ggplot(aes(scrublet_call, nCount_RNA)) +
  geom_boxplot() +
  stat_compare_means()

scrb.res |>
  ggplot(aes(scrublet_call, nFeature_RNA)) +
  geom_boxplot() +
  stat_compare_means()

# worker number = sample number: ~2.5 min with 11% doublet
# default MulticoreParam() ~ 2.5 min with 11.2% doublet
# SerialParam() ~6 min with 11.2%
# SnowParam() ~4 min with 10.9% & warning
system.time(
  res.dbl <- sobj_liv |>
    as.SingleCellExperiment() |>
    scDblFinder(clusters = 'seurat_clusters', samples = 'orig.ident',
                returnType = 'table', BPPARAM = SnowParam(log = T)) |>
    as_tibble(rownames = '.cell') |>
    filter(src == 'real') |>
    mutate(.cell = .cell, scdbl.class = class, .keep = 'none'))

res.dbl |>
  summarise(dbl.rate = sum(scdbl.class == 'doublet') / n())

# encapsulated SerialParam: ~6.7 min with 10.9%
# encapsulated MulticoreParam with log=T: ~3 min with 10.6%
# encapsulated MulticoreParam with log=T & workers=6: ~2.5 min with 10.8%
system.time(
  res.dbl <- sobj_liv |> mark_doublets2())

sobj_liv %<>%
  left_join(res.dbl)

sobj_liv |>
  ggplot(aes(seurat_clusters, fill = scdbl.class)) +
  geom_bar(position = 'fill')

sobj_liv |>
  ggplot(aes(scdbl.class, nCount_RNA)) +
  geom_boxplot() +
  stat_compare_means()

sobj_liv |>
  ggplot(aes(scdbl.class, nFeature_RNA)) +
  geom_boxplot() +
  stat_compare_means()

# doublet hard threshold in Zhang 2020
# ~3.9%
sobj_liv |>
  summarise(dbl.rate = sum(nCount_RNA > 25e3 | nFeature_RNA > 4e3) / n())

# doublet hard threshold in Zhang 2021
# ~1.1%
sobj_liv |>
  summarise(dbl.rate = sum(nCount_RNA > 4e4 | nFeature_RNA > 5e3) / n())

sobj_liv |>
  filter(scdbl.class == 'singlet') |>
  write_rds('CRC-I/data/li_2021_liver/li2021crc.rds')
